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Abstract 

The density of states near zero energy in a graphene due to strong point defects with random 
positions are computed. Instead of focusing on density of states directly, we analyze eigenfunctions 
of inverse T-matrix in the unitary limit. Based on numerical simulations, we find that the squared 
magnitudes of eigenfunctions for the inverse T-matrix show random-walk behavior on defect posi- 
tions. As a result, squared magnitudes of eigenfunctions have equal a priori probabilities, which 
further implies that the density of states is characterized by the well-known Thomas-Porter type 
distribution. The numerical findings of Thomas-Porter type distribution is further derived in the 
saddle-point limit of the corresponding replica field theory of inverse T-matrix. Furthermore, the 
influences of the Thomas-Porter distribution on magnetic and transport properties of a graphene, 
due to its divergence near zero energy, are also examined. 

PACS numbers: 81.05.ue, 61.72.J-,71.15.-m 
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I. INTRODUCTION 



Recently, the isolation of single-layer graphene^ has revived much interest in studying 
two-dimensional (2D) Dirac fermions. One of the peculiar properties associated with 2D 
Dirac fermions is the unusual electronic properties in the presence of defects and disorders. 
In the context of cuprate superconductors, where quasi-particles are also 2D Dirac fermions, 
disorders have masked the d-wave nature and hindered its discovery. It was later realized 
that point defects may change the density of states (DOS) near the Dirac point and strong 
point defects may even induce quasi-localized states or magnetic moments near zero energy 
in d-wave superconductors^. In the case of graphene, it is found that there is finite density 
of states due to weak disorders^. For strong disorders, it was observed that ferromagnetic 
state can be induced by bombarding a graphite with protons^. The induced magnetism is 
further confirmed to be resulted from vr-electrons^. This fact, together with recent obser- 
vation of ferromagnetism in disordered graphene^i^, shows that graphene with defects could 
become ferromagnetic. In addition to magnetism, graphene also reveals anomalous transport 
properties in the presence of strong disorders, where in the presence of vacancies, instead 
of decreasing, the conductivity is found to increase^. These observations clearly indicates 
that in the presence of strong disorders, 2D Dirac fermions may behave very differently from 
what is expected for clean or weak disordered graphene. 

Experimentally, there are many possible forms of disorders in graphene^. For large de- 
fects such as cracks, they tend to contain the so-called zig-zag edges, where localized states 
would appear near the edge^^ and induce magnetic behavior—. In this case, magnetic mo- 
ments arise from localized states and interact via RKKY (Ruderman-Kittel-Kasuya-Yosida) 
interaction, which tends to make graphene antiferromagnetioi^. Hence the most possible 
candidates for the observed ferromagnetism in graphene are defects of small sizes or sim- 
ply point defects. Here the simplest point defects are single-atom vacancies or hydrogen 
chemisorption defects. These kinds of defects generally create complicated disturbances in 
graphene and may even form ordered structures^. However, for low density of quenched 
defects, they can be simulated by a large potential m on a lattice point without distortion 
of nearby lattice points^. 

Theoretically, extensive studies on a single defect have been performed on d-wave 
superconductors^. It is known that a zero-energy electronic state would arise near a point 
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defect with m — )■ oo or a circular disk in a 2D Dirac Hamiltoniani^. Furthermore, the elec- 
tronic wavefunction is semi- localized with amplitude decaying as 1/r at large distance r—i^^. 
The semi-localized behavior is clearly revealed in the observed STM images of long-range 
(v^ X v^)-R30° superstructure in grapheno^i^. For finite density of defects, one expects that 
semi-localized electrons interact strongly and may form an impurity band^»^. Nonetheless, 
conflicting results based on either perturbative or non-perturbative approaches are reported^. 
The residual DOS near zero energy is predicted to be either finite^^, infinite^ii^, or vanishing 
with different power laws in energy^. This issue remains unsolved. 

While quasi-particles in both cuprate superconductors and graphenes are 2D Dirac 
fermions, the situation is quite different for graphene. For neutral graphene, even though 
excitonic effects are expected to be large^, for low energies and large distances, the screened 
Coulomb interaction is shown be long-ranged^i with renormalized dielectric constant. Fur- 
thermore, the electron itself is the quasi-particle and carries a definite charge. These differ- 
ences make graphene behave totally different from that of cuprate superconductors in the 
strong disorders. In particular, without being masked by superconductivity, direct manifes- 
tation of the impurity band is possible in graphene. Therefore, investigation on graphene 
with strong defects would provide an unique opportunity to clarify the issue of DOS near 
zero energy for 2D Dirac fermions with strong disorders. This is recently pointed out in 
Ref . [igl . In that paper, the wavefunction for finite density of defects is constructed. By us- 
ing the wavefunction for two defects, it is shown that ferromagnetic state is favored for large 
distances between two defects. However, for finite density of defects, the problem of finding 
DOS is mapped to an equivalent problem of finding the DOS of a random matrix. One has 
to assume that the matrix elements are independent random numbers to demonstrate the 
induced ferromagnetism?i2. While the predicted DOS (Wigner semi-circle law) appears to be 
consistent with results obtained by self-consistent Born approximation^^, to confirm that the 
observed ferromagnetism and anomalous transport properties of graphene are consequences 
of the impurity band, one needs to go beyond self-consistent Born approximation and to 
resolve the issue of how the DOS of 2D Dirac fermions changes in strong disorder limit. 

In this paper, we re-examine the density of states of a graphene due to strong point 
defects. In particular, we show that the inverse T-matrix for Nj point defects can be exactly 
mapped to a. Nj x Nj symmetric Euclidean Random Matrix in which one cannot treat 
the matrix elements as independent random numbers. Instead of focusing on the DOS 
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directly, we analyze magnitudes distribution for eigenfunctions for the derived Random 
Matrix. Remarkably, we find that squared magnitudes of eigenfunctions show random- 
walk behaviors on defect positions. As a result, the distribution of squared magnitudes of 
eigenfunctions for the Euclidean Random Matrix follows the Porter-Thomas distribution. 
Further analysis shows that eigenvalues (A) of the corresponding Euclidean Random Matrix 
also follow the Thomas-Porter distribution^^ and the DOS near zero energy for infinite u is 



Here n/ is the density of defects and (|A|) is the average of |A| over defect configurations. 
X{E) is given by \{E) = —^^E In \E/D\ with D = 3t/2 and t being the hopping amplitude 
of the electron. This form of the density of states is valid when \E\ t and we found that 
(|A|) ~ y/nj shows random- walk behavior. The resulting density of states has strong effects 
on magnetic and transport properties of graphene. We re-examine the effect of the long- 
range Coulomb interaction with renormalized dielectric constant and show that the resulted 
DOS supports ferromagnetism for any finite density of defects. At finite temperature, the 
linear extrapolation of magnetization curve indicates that ~ 600 — 700-ft", in agreement 
with experimental observations. 

This paper is organized as follows. In Sec. II, we lay down the theoretical formulation and 
show that the inverse T-matrix for Nj point defects can be exactly mapped to a Nj x Nj 
Euclidean Random Matrix. In Sec. Ill, we use both analytic arguments and numerical 
simulations to derive the density of resonant states. In Sec. IV, we reexamine effects of 
the screened long-range Coulomb interaction. We show that the competition between the 
exchange energy and kinetic resonant energy leads to ferromagnetism for infinite on-site 
potentials. The magnetizations both at zero and finite temperatures are also calculated. In 
Sec. V, we conclude and discuss possible effects for weak impurities. Appendix A is devoted 
to more rigorous derivation of the Porter-Thomas distribution in the saddle-point limit. 

II. THEORETICAL FORMULATION 

We start by setting up the framework for investigating the effects of defect. It is known 
that electrons in the vr band of an infinite graphene can be well described by a tight-binding 
Hamiltonian i^cr. As shown in Fig. [T|, the lattice of graphene is bi-partite. If we label 




(1) 
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the bi-partite lattice points by A and B, Hq consists of hopping only for nearest A and B 
with a hopping amplitude t. Hence if defects are located at fi with i — 1,2,3, - ■ ■ ,Ni, the 
wavefunction for an electron then satisfies 



^ Ni \ 

K i=l J 



(2) 



Here and in the following, both r and fi are restricted to points on the honeycomb lattice 
shown in Fig. 1. To find the effects of defects on the electronic state, it is sufficient to 




FIG. 1: Honeycomb lattice of graphene. The lattice is bi-partite, labeled by A and B , with hopping 
amplitude between nearest A and B being t ~2.7eV. The lattice constant a = 2.46A is the distance 
between two nearest B points. 



calculate the Green's function G{f, f , E), which describes the amplitude for the electron to 
propagate from r* to r and satisfies 



(3) 



where H = Hq + u Ylf=i ^f,fi- For clean graphene, the Green's function will be denoted by 
G°(r, f, E). In the Fourier k space, it is convenient to reorganize the wavefunction into if) a 
and ■0B for A and B sublattices. Then G^{k) is the inverse of the 2x2 matrix, E-\-iQ'^ —HQ{k), 
with Hglk) being given by 

A{k) 
A*{k) 



Hoik) 



(4) 
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where A(fc) = —t[2e''''y°-'^^ cos{kxa/2) + e *'=!/'^/v^]. More explicitly, one finds 
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1 1 



E + ^0+ - |A(A;)| E + iO+ + \A{k)\ 



(5) 



(6) 



In real space, it is more convenient to use lattice vectors r to carry indices for A and B 
sublattice. Therefore, G is no longer a 2 x 2 matrix and can be expressed in terms of as 

Ni 



G{f, r", E) = G°(r, 7^,E) + uJ2 f,, E)G{fi, r", E). 



(7) 



i=l 



Clearly, to find G{f,r',E), one needs to find G{fi,r',E) in Eq. ([7]). For this purpose, one 
sets r to fi with i = 1, 2, 3, ■ ■ ■ ,Nj in Eq. ([71) and solves (^(ri, r*, E) in terms of G''(ri, r'). 
If we replace the notation Gifi^r'^E) by Gr^.r* with E being suppressed, we obtain 



l/u 



\ 



-Gl 



rNjA 



GO., . l/u - G% 



llu-G% . / 



I Gi „ \ 



Here the matrix on the right hand side is the T-matrix whose inverse determines resonant 
energies and can be separated into real and imaginary parts 
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y-1 



( l/u - Gn -Qi2 

—Q21 i/u — G22 

—Q3I ~Q32 



-QlNi 

-Q2NI 
-Q3NI 



( hi /12 

hi h2 
hi h2 



hNr ^ 

hN, 



(9) 



l/u — QnjNi j \ IniI Ini2 ■ ■ IniNi J 



where Qy and hi are the real and imaginary parts of G^-. Note that due to the Kramers- 



Kronig relation, Qy and hj are related by 



(10) 



Therefore, real (T^^) and imaginary parts (T^ ^) of T ^ can be diagonalized simultaneously. 
In particular, their eigenvalues are also related by the Kramers-Kronig relation 



(11) 



It is thus clear that the resonant energies of the Green's function G are determined by zeros 
of eigenvalues of Tr. Therefore, resonant energies due to defects are determined by 



1/u-Qn -G12 
-Q21 l/u -G 22 

— Q3I ~Q32 



-QlN, 

-Q2N, 
-Q3NI 



l/u - Qni 



Ni 



0. 



(12) 



Note that the above condition is exactly the same as the one obtained via the constructed 
wavefunction for defects^ and should be compared to the similar equation obtained in the 
context of d-wave superconductors^-. Since a graphene without defect is translationally 
invariant, one has G^^ = G^{fi — rj). Therefore, diagonal terms in Eq.f ll2p are identical and 
are equal to \{E) = l/u — ReG^{0,E). Hence if the positions of defects are random, solving 
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Eq.( !T2|) is equivalent to finding eigenvalues of the random matrix 



Hr 



( ° 


Q12 ■ 


■ QlNi ^ 


G21 
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Q31 


Q32 ■ 


■ Q3N1 



(13) 



\ Qnii Qni2 ■ ■ y 

Furthermore, if one defines the density of eigenvalues for Hj by 



P(A) 



(14) 



with M being the total number of lattice points and A„ being the n-th eigenvalue, the density 
of resonant states is given by 

d\{E) 



D{E) = V{X{E)) 



dE 



(15) 



Therefore, it is sufficient to find the distribution of eigenvalues for Hj. We note in passing 
that if values of Xq form a band after averaging over defect configurations, it implies that the 
averaged (Ac) is independent of E. Eq.l fTTl) then implies that except for contributions from 
diagonal terms off-diagonal terms do not contribute to the imaginary part of eigenvalues. 
Hence if values of A form a band, one has T^^ = — ImG"'(0, E)I. Since ImG^(0, E) oc E, this 
result implies that the inverse of lifetime for resonant states is proportional to E, consistent 
with experimental observational. 



III. DENSITY OF RESONANT STATES 



In the last section, it is shown that the density of resonant states is determined by the 
spectrum of Hj. Since each element, Qy, depends on positions of defects, they fluctuate 
randomly. In the simplest approximation, one treats each element as an independent ran- 
dom number. The density of states is characterized by the Wigner semi-circle law^S. As 
indicated earlier, this approximation appears to be equivalent to the self-consistent Born 
approximatioii^^. A closer examination of Hj shows that the dependence of each matrix 
element on the position makes them correlated. Hence one cannot treat each element as 
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an independent random number. Indeed, it was realized in different context by Mezard et 
al.— that such random matrices form distinct classes known as Euclidean Random Matrices, 
whose spectrum depends on the functional form of the matrix element on rj . 

It is generally difficult to find the exact spectrum for any given Euclidean Random Matri- 
ces. For defects on graphene, however, it turns out that the spectrum of Hi follows a simple 
form known as the Porter-Thomas distribution^^. In this section, we shall focus on the study 
of the spectrum by numerical simulation. An analytical derivation based on saddle-point 
approximation will be relayed to the Appendix. 

We start by noting that since one expects that the energies of resonant states are close 
to zero, as a first step, we can approximate each matrix element by Qij{E = 0). We shall 
see later that the error due to this approximation is small for ~ 0. In this approximation, 
by using Eqs.(l5]) and ([6]), one finds QAA{,f,E = 0) = QBsi^.E = 0) = and GsAU^i) = 
G\^{i,j). Hence Hj is a symmetric matrix. Furthermore, since in the second quantization 
form. Hi = J^ij^ij'^Ai'^Bj + h.c, we find that Hj goes to —Hj under the particle-hole 
transformation: c\- — )■ — c^j and cbj — )■ c^bj- Therefore, the spectrum is particle-hole 
symmetric, i.e., P(— A) = 'D{X). In addition of being particle-hole symmetric, Hj itself 
also supports energy states exactly at zero energy due to the unbalance in the number of 
lattice points in A and B^^. Since the number of zero energy states is equal to \Na — Nb\, 
if lattice points are randomly assigned to A or B, one finds |A^^ — Nb\ ~ and hence 
their contribution is negligible in the limit of M — )■ oo with Nj/M being fixed at the defect 
density nj. Therefore, in the following, we shall focus on density of resonant states for the 
case with Na = Nb to avoid complications due to extra zero energy states. 

For high density of defects, because the positions of defects sample sufficient lattice points, 
Hj can be diagonalized by Fourier transformation. Hence eigenvalues of Hj are proportional 
to the Fourier transformation of Qy. We find that 

V{\) oc /" /" ^ Isfx- TT^) +5(\ + r-^] . 

^ ^ J J (2vr)2 [ V |A(q)i; V |A(q)lA 

(16) 

In this case, because < |A(g)| < 3t, we obtain A > l/3t. Therefore, there is no resonant 
defect state near zero energy for sufficient high density of defects. 

For low density of defects, the separation between any two defects is large. In this case. 
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FIG. 2: (Color on-line) Participation number {p = 1/ ^'^d histogram of |</>aP of different 

A's and p's (indicated by sub-indices of A) for a typical defect configuration simulated with Nj = 
1000 and M = 1000 x 1000, i.e., n/ = 0.001. Here red sohd hues are the fitted Boltzmann 
distributions. 
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by using Eqs.(l5]) and ([6]), we find that for < \E\r/v ^ 1- 



GaaIt, E) = GsBir, E) = cos ( — 1 Eln ■ 



2ttv^ 



3a 



QAB[r, E) = QbaKt, E) = sm 



2ttv r/a \ 3a 

Here v = 3ta/2. While for = 0, as we indicated earher, QAA^r, E = 0) = GsBir, E = 0) 
but QAB^r, = 0) is given by Eq. (|T8l) . For r = 0, we obtain 

^2 



(17) 
(18) 



gAA{o,E) = gBB{o,E) 



v/3a^„, a\E\ 



-E\n 



(19) 



To motivate it, instead of focusing on DOS directly as done in the d-wave 
superconductors^'^^, we analyze the distribution of the eigenfunction amplitudes (pxifi) of 
Hi at a fixed eigenvalue A 



Here 4'\{ri) is normalized so that 



(20) 



1. 



(21) 
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Hence if tliere is no bias on partitioning one expects |0p follows the Boltzmann type 
distribution, P(|0p) oc e""'''^'^. Indeed, in the limit M oo. Porter and Thomas^ derived 
the following distribution 



where s = |</)p and (s) is the average of The same distribution can also be derived in 
the non-linear sigma model-^. The Porter-Thomas distribution, however, is not universal 
and is valid only when the system is sufficiently chaotic^-. Since the matrix element Qij 
decays slowly (1/t), 0a at each point fi is determined by all other defects with random 
positions. In other words, Hj is a random hopping model in which 0j characterizes density 
of random walkers on defect position r,. Since the probability for finding a random walker at 
the traveling distance r is proportional to e~''^/^^'"^\ by comparison with Eg. (1221) . one expects 
that the Porter-Thomas distribution works for Hj with playing the role of distance. More 
explicitly, for a random walker described by r{t), one finds (r^) oc t at time t. Here t 
characterizes the number of attempts in a random walk. By analogy, Nj would be the 
number of attempts. Therefore, we expect 



Based on Eq. (fT8|) . we perform extensive numerical analysis on the statistics of eigenstates 
of Hj. To see if there is correlation between distribution and localization of 0a, we also 
analyze the participation number p = I0(^i)l^ ^^e distribution for different 

participation numbers. Fig. |2] shows the statistics of wavefunction amplitudes for a typical 
defect configuration. It is clear that regardless of whether the eigenfunction is localized or 
not, distribution of amplitudes follow the Porter- Thomas distribution for all participation 
numbers. 

For different A, in addition to Eq.( l2T|) . partition of eigenfunction amplitudes has an 
addition constraint 




(22) 




(23) 




(24) 




distribution for A also follows the Porter-Thomas distribution 




(25) 
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Here according to Eq. (l23|) . we expect (|A|) oc ^Jnl. The proportional constant will be deter- 
mined numerically. The normalization T> in Eq. (l25|) is chosen by requiring dXD{X) = rii. 
Note that for later use in the calculation of magnetization, the normalization of V has to 
be done by taking into account the presence of Dirac band. 

Fig|3](a) shows a typical spectrum of our numerical simulations of the spectrum averaged 
over 1000 defect config well described by the 




FIG. 3: (Color on-line) (a) Averaged spectrum of Hj (with Qij{E = 0) as the matrix element) 
over 1000 defect configurations. Here n/ = 1000, M = 1000 x 1000 and we have set t = 1. 
Black circles are numerical results while the red line is the fitted Porter-Thomas distribution with 
P(A) = 2.17e~'^^'^^^l'^l /y^|A[. Inset: The corresponding density of electronic states for u = oo. 
(b) Random-walk behavior of cj): The dependence of (A) on the defect density y/nj shows linear 
behavior with a slope 0.464. Here open circles are numerical results obtained by the fitted Porter- 
Thomas distribution while the red line is the linear curve of slope 0.464. There is a small error 
offset by -0.0038. 

Porter-Thomas distribution. Fig|3]^b) shows the fitted parameter (A) versus density of de- 
fects. It indicates that (A) follows a simple form of ni by (A) ~ ^JnJ- Once one knows the 
spectrum of Hi, by using Eq.(!T9ll. the density of resonant energies can be found by setting 
A = \/u — Qaa/bb{S^, E). This results in Eq.([I]). In the inset of Figj3]^a), we show the 
corresponding electronic DOS for u = oo. It is clear that the DOS diverges at E = 0. 

We close this section by checking the validity of setting E = in Qij{E). For a given 
finite E, because A = 1/u — Q{0, E), there is only one value of A corresponding to the given 
E. Hence for a given E, only the spectrum at A = 1/u — Q{0,E) is correct. To get the 
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whole spectrum, it is necessary to vary E and obtain the spectrum at each X{E) one by 
one. Note that by using Eqs.([5]) and IQ, one finds that Qij{—E) = —Qij{E) and hence the 
resulting spectrum is still particle-hole symmetric. In FigJH we show the comparison of the 
spectrum for Hj by using Qij{E) and Qij{E = 0). It is clear that the difference is small and 
both spectra follow the Porter-Thomas distribution, in agreement with the derivation in the 
Appendix that is based 

160 
120 

Q 



" 0.002 -0.001 0.001 0.002 

X 

FIG. 4: (Color on-line) Comparison of the spectrum of Hj determined by using Qij{E) (open 
circles) and Qij{E = 0) (red crosses) as matrix elements. Here n/ = 1000, M = 1000 x 1000 and 
we have set t = 1. Both spectra can be fitted with the Porter-Thomas distributions with slightly 
different (|A|). 

IV. EFFECTS OF COULOMB INTERACTION AND FERROMAGNETISM 

In this section, we discuss effects of the Coulomb interaction due to the change of density 
of states. To include the effects of Coulomb interaction, we note that for neutral graphene, 
even though excitonic effects are expected to be large^, for low energies and long distances, 
screening can be taken into consideration by the renormalization of v and the dielectric 
constant e^^. Hence for low density of defects in which separation between any two defects 
is large, one needs to replace v by vr in Eqs. ( 1T71) and ( ITSi) . This would effectively replace 
the hopping amplitude t by t^. 

As indicated in the introduction, strong disorders in a graphene are a possible source 
for the observed ferromagnetism. To examine whether the Porter-Thomas type distribution 
supports ferromagnetism, we first note that the normalization adopted in Eg. (125!) has to 
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be corrected by taking into account the conservation of states. As indicated in Figj5l since 
resonant states replaces states in Dirac band, it requires a cutoff A in the impurity band 
so that numbers of states for the impurity band and the Dirac band are equal. Since 

D(E) 




-A A 



FIG. 5: Schematic plot of the impurity band and the original Dirac band. Here the solid lines 
at center represent the impurity band when u = oo. For a small finite u, the impurity band 
(represented by the dash line) is shifted into the Dirac band and disappears 

number of states for the impurity band per site is nj, integration of the Dirac band yields 
A = ^/Tr/Avnni. By including the cutoff A, appropriate normalized V\ is given by 

rij [AL 



-e 2(|A|) _ 



(26) 



/M^erf(v/AA/2(|A|)) 
Here erf is the error function and Aa = A (A). Note that when v is renormalized to vr, both 
(|A|) and A a are renormalized by the same factor vr/v. 

To investigate the magnetism, we note that the electron wavefunction is related to 
the eigenfunction (p of Hj as follows^ 

Mr^)= J2 ^^^^E- (27) 

j =defect positions 

Here A'^ is proportional to (f)j. The normalization of A'^ is determined by (Xli^el^i)) ~ 
1. The applicability of the Porter-Thomas distribution implies that ipj (thus A'^) follows 
Gaussian statistical^. Hence we have 

{A%Aj,) = r6.,. (28) 

By expressing Qij = jij2^q'0{.(l)^^^^^''~^^'' and using the fact that QAA{r,E = 0) = 
gBB{r, E = 0) = 0, we find F = l/(A^/7) with 



7 = ^AB{q)GAB{-q)- 



(29) 
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We shall include the Coulomb interaction by calculating the exchange energy. For a neutral 
graphene, it is known that screened Coulomb interaction is still long-rangeciii2^ 

= 1^ E CtC„^C],C,,., (30) 

where e is the renormalized dielectric constant and is roughly 2.3eo. To obtain the exchange 
energy, Eq.(l27]) is replaced by Cj^ = Y^ej '^EGiji^)'^Ea- By setting any pair of CLCi^v' 
by its average value {Cl^CE'a'), using the fact (A^^^Aj^j^Al^A^^) = (A*^^A^J(v4^^A|J and 
approximating Qij{E) by Qij{0), we find that the exchange energy is given by 

Eex = , B 31) 

with ^ 

B = Y.T^\{fl^-^^^^^ ' (32) 

where = N^r/Nj are fractions of electrons in the spin state a. By approximating Qij{E) 
by Gij{0) and expressing Qij in Fourier space, we find 

B _ IGnnj 

^. GAB{q)GAB{-q)GAB{q)GAB{-q) (33) 

q,q' ^ 

where only the i and j in the same sublattice would contribute. Since for E' ~ 0, Eq.(l6]) 
implies Gab diverges near Dirac points qu. The main contribution in the integral of B comes 
from regions of qo- By setting q to any one of the Dirac points in the factor l/\q + q'\, we 
find 

-nh'. (34) 



M 4a 

In the ferromagnetic state, we have ^ with E^r being the corresponding Fermi energy 
for the spin state a. The net spin moment is proportional to m = n-|- — rii. Substituting 
Eq.f p^ back to E^x, we find that the exchange energy per site due to m is given by 

E,x 6^2 + ^3 2 /^KN 
= Tij. (35 

For an undoped graphene, E^ = —E^. In this case, the net spin m can be expressed as 
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FIG. 6: (Color on-line) Magnetization per defect for u = oo with screened (e = 2.3eo, = 1.3v) 
and unscreened (e = eo, f/j = v) Coulomb interactions, (a) Magnetization versus defect density at 
zero temperature (b) A typical temperature dependence of magnetization for ni = 0(10^^). There 
is no sharp transition temperature. However, by linear extrapolation, one finds that Tc is around 
600 - imK. 

m = 2 j^^ dED\{E), while the change of the total energy in the impurity band per site is 
A/c = 2 j^^ dEED\{E). The minimization of AA; + E(,x/M with respect to E'^ then leads to 



(2 + V3)e 



IGvrea 



/ Vj,{\)dX = E^. (36) 
Jo 



Solving Eq. (l36l) yields E^ which in turn determines the magnetization per defect at zero 
temperature. The same calculation can be easily generalized to any finite temperature T. 
In this case, the magnetization is still determined by the minimization of Ak + Eex/M with 
respect to E"^ except that now m and Ak are replaced by 

/oo 
dEDK{E)[n^{E)-n^{E)], (37) 
-oo 

/oo 
dED^{E) [n^{E) + n^{E) - 2no{E)] , (38) 
■oo 

where n^^E) = 1/ (^e^(^+^'^^ + 1) are the Fermi-Dirac distributions for a =t or J, and no{E) = 
l/(e^^ + 1) with /3 = 1/kBT. 

In Figjnt^a), we show the magnetization at zero temperature for u = oo with screened and 
unscreened Coulomb interaction by solving Eq. (l36|) . It is seen that screening reduces the 
magnetization. Furthermore, due to the divergent DOS at E = 0, ferromagnetism persists 
down to zero defect density and magnetization increases as defect density increases. FigJHI^b) 
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shows a typical temperature dependence of magnetization for rij = O(10~^). The temper- 
ature dependence shows a quasi-hnear behavior with a Boltzman tail. To compare with 
experiments, we perform the linear extrapolation of magnetization curve, which indicates 
that Tc ~ 600 — 700i^, in agreement with experimental observations^. 

V. SUMMARY AND DISCUSSION 

In summary, in this work we have shown that in the strong disorder limit, a resonant 
impurity band is induced in a graphene. By combining analytic arguments and numerical 
calculations, we show that the density of resonant states is governed by the principle of 
equal a priori probabilities for squared magnitudes of eigenfunctions of a Euclidean Random 
Matrix. For large on-site defect potential, the principle of equal a priori probabilities shows 
that the density of resonant states is characterized by the Thomas-Porter distribution and 
is divergent near zero energy. Furthermore, we show that the observed ferromagnetism 
is due to the combination of strong disorder and long-range Coulomb interaction. The 
linear extrapolation of magnetization curve indicates that Tc ~ 600 — 700K, as observed in 
experiments. 

In addition to the magnetism, the impurity band enhances the transport^^. This is 
consistent with experimental observations^ but is quite different from ordinary impurity 
states even though in the calculated participation number of Hj in FigJ2l some eigenfunctions 
are localized. The crucial difference lies in the semi-localized nature of the electronic states 
as revealed in Eq. fl27|) . Here even though A-'^ (thus (p^) is localized, due to that Q ~ 1/r, 
ipE will not be exponentially localized around defect positions. The participation number 
for tpE itself is of the order of (InM)^, indicating its semi-localized nature. 

While so far in this work we only consider the strong disorder limit, the results also pro- 
vide some insight into the weak disorder region. As illustrated in Figl^l for weak disorders, 
u is small, the impurity band is shifted into the Dirac band. In this case, while the majority 
weight of the impurity band disappears, its tail still sweeps through zero energy and con- 
tributes small but finite DOS. As indicated above, these density of states generally enhances 
the transport. This explains why when graphene is made cleaner, the conductivity, instead 
of increasing, decreases and appears to approach an universal constant^. While the impurity 
band cannot account for the exact value of the universal conductivity, our results serve as a 
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useful starting point for obtaining corrections to the conductivity. 
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Appendix A: Saddle-point limit and Porter-Thomas distribution 

In this appendix, we shall show that the eigenvalue distribution for Hi follows the Porter- 
Thomas distribution in the limit oi Ni ^ oo but with the defect density rti being fixed. We 
start by noting that the spectrum of Hi can be found by calculating the resolvent 

where M is the total number of lattice points and (■) is the average over the random config- 
urations of defects. Clearly, we have ^^(A) = — ^lm/?(A + iO^). As shown in the text, since 
the spectrum has particle-hole symmetry, we shall set A to |A| and consider only positive 
A. The evaluation of R{z) can be reformulated by a replica field theorj*^^ via the following 
identity 

Riz) = lim^^(e-"^'^i°s(^-^^)) 

^ hm^^/ / \. (A2) 

The term 1/ det(2; — Hif can be re-expressed by n rephca complex fields 0a (a=l,2,3,..,n) 
as follows 

Ni n 



1 



det(;z - HiY 



i=l a=l 



(A3) 



Up to now (pa is only defined on defect sites. To remove this constraint, one introduces the 
field ipa defined on every lattice site and impose 6[ipa{r) — ^"(^r.rj- The constraint of 
the delta function can be removed by using the identity 5{F) = f dipa^''^"'^ ■ Here ipa is the 
replica field. After integrating out (f)°; and ipa, the resolvent can be expressed as2^ 

^W = -lim4.^, (A4) 
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where the partition function Z is given hj Z = J Vipe ^ with S being given by 

n 

S = -J2Y1 ra{n)F,ji^a{r^) -niY, e~-^ (A5) 

a=l i,j i 

Here Fij = Q^^ — 5ijQl^^ . Note that because = 1/{E — H), we have that for ~ 0, 
Fij = —{HQ)ij. In other words, —F has the same form as the tight-binding Hamihonian for 
graphene except that it only acts on defect sites. 
After substituting Z back to i?, we find 

K[z) = — hm ■ t: 

^ ^ n^o nM f Vipe-^ 



hm 



M / Vtpe- 



s 



(A6) 



Here we have made use of the equivalence among different replica component a and the 
equivalence among different positions in the second equality. 

It is clear that in Eq. flA6p . the replica symmetry is broken. One needs to perform inte- 
grations for ipi and ipa with a 7^ 1 separately. For ipi, the integrand can be rewritten as 
J2i with SI given by 

Sl=^-^l^^-\n\Mm' + S. (AT) 

Since we shall be interested in z ~ 0, i.e., energy near zero, in the saddle-point approxi- 
mation, integration over ipi is dominated by the maximum of SI, which is determined by 
■rrr^T^Sl = for all rt-. We find that maximum of 5? satisfies 

z 

-V^?*(rl)5^F,,AV.) = 0. (A8) 
j 

It is clear that for low density, we can expand ip^ in term of rij. We find \tpf{fi)\'^ / z = 

l—nje~^^^ ^'^"^''^^^^ — ip^* (fi) J2j Fijipiifj) H . Because Fij is finite, we obtain ipiifi) ~ ^/z 

for all Tj. As a result, the integration of ipi in Eq. ( ]A6I) can be approximated as 



i[ri e - 
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Here 6 is a constant that results from Fij and (SI)" is the 2nd derivative of 5*^ with respect 
to ipi and is proportional to 1/z. li is the integration over 'ipaifi) with a ^ 1 and is given by 

= (e-iJ^^^il'^»(^^')l')'Z;_i, (AlO) 
where ()' is the average with respect to 5" and Z'^^_-^ = J Vipa^'^' with 5" being given by 

n 

S' = -5^5^^:(r-;)F,,^,(r-;-) -n,e-i5^e-^^''l'^"(^^-)l'. (All) 

a=2 i,j i 

It is clear that different a and i are equivalent. Therefore, we obtain 

^$:/. = (e-^l^^(-)lV^U. (A12) 

i 

Combing Eqs. (IA9p and flA12p gives the limiting behavior of the numerator for small z. To 
obtain the spectrum for small A, one needs to find the analytical continuation by replacing 
2; by A + iO^. Clearly, the factor e~^^/-\/i only contributes the real part. Together with the 
fact that the denominator is Z = ^-^'^^^"si^-Hi) ^ which goes to one when n approaches zero, 
we find that 

V{X)=nj^- lim Im(e^^^)'Zii. (A13) 

A-s>0,n-s>0 

Both (e^''^'''-''*-''^)' and Z'_^ can be calculated perturbatively with finite results in the limit 
A — i- 0^^. After appropriate normalization, one finds the spectrum of A follows the form 
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